%Structual estimation (main file), including decreasing sunk costs



clear
tic

choice=1:13;
beta=0.95;
T=100;
length=0:T;
df=beta.^length;
t0=2006;
tN=2009;
tbar=tN-t0+1;
%AI_fs=[51:63,65:66,68:72,76:94,96:119,121:145,147:150];


start=1;
page=start:(start+29);
narray=numel(page);
%narray=s;
%coef=zeros(narray,2);
%ste=zeros(narray,2);
%fv=zeros(narray,1);

fs0=dlmread('fs1.csv');
fs0=sortrows(fs0,[1,2,5]);
ddom0=dlmread('mks1.csv');
ddom0=sortrows(ddom0,[1,2,3]);
FSTA=zeros(size(fs0,1),size(fs0,2),narray);
DOM=zeros(size(ddom0,1),size(ddom0,2),narray);

% Construct other inputs from "regvar*"
bd=fs0(:,3,1);
%nobd=bd(:,ones(T+1,1));
%bdtype=(bd<=25&bd>0)+(bd<=80&bd>25)*2+(bd<=193&bd>80)*3+(bd>193)*4; %bdtot: 25%: 25; 50%: 80; 75%: 193
tmpbd=unique(bd);
tmpbd(tmpbd==0)=[];
bed25=prctile(tmpbd,25);
bed50=median(tmpbd);
bed75=prctile(tmpbd,75);
bdtype=(bd<=bed25&bd>0)+(bd<=bed50&bd>bed25)*2+(bd<=bed75&bd>bed50)*3+(bd>bed75)*4;
%bdtype=[(bd<=bed25&bd>0), (bd<=bed50&bd>bed25), (bd<=bed75&bd>bed50), (bd>bed75)];



%### UPDATED on 05/07/2020:Calculate the mktdom at the beginning here NOT using the one in the policy
%function because we need the endogenous mktdom (including own choice) so
%each hospital in the same market has the same mktdom
textFileName26=sprintf('regvar6_sing.raw');  %the input file  
textFileName27=sprintf('regvar7_sing.raw');  %the input file     
textFileName28=sprintf('regvar8_sing.raw');  %the input file     
textFileName29=sprintf('regvar9_sing.raw');  %the input file    

char6=dlmread(textFileName26); %the matrix of characteristics of stand-alone hospitals: ahaid year bdtot prechs adj_vnd hsa
char7=dlmread(textFileName27);    
char8=dlmread(textFileName28);
char9=dlmread(textFileName29);

allchar=[char6;char7;char8;char9];
allchar=sortrows(allchar,[1,2]); 
pre=allchar(:,4);
act=allchar(:,5);  %extract the current action to compute the probability of the actual choice


%{
nowDOM=zeros(size(allchar,1),numel(choice));
for i=1:size(allchar,1)
    tmpmkt=allchar(allchar(:,end)==allchar(allchar(:,1)==allchar(i,1)&allchar(:,2)==allchar(i,2),end),4);
    tmpmkt(tmpmkt==1)=[];
    tmpmkt(tmpmkt==2)=[];
    tmpmkt(tmpmkt==13)=[];
    if isempty(tmpmkt)==0
       occur=histc(tmpmkt,choice)';
       occur=reshape(occur,1,[]);           
       mktdom=choice(occur==max(occur));
       nowDOM(i,mktdom)=1;
    end    
end
nowDOM=[allchar(:,1:2),nowDOM];
%}


% Define cuboid for fminsearch
df_cub=repmat(df,size(FSTA,1),1,narray);
bd_cub=repmat(bd,1,T+1,narray);
bdtype_cub=repmat(bdtype,1,T+1,narray);
gg_cub=zeros(size(DOM,1),T+1,narray);
fchs_cub=zeros(size(FSTA,1),T+1,narray);
ind_cub=zeros(size(FSTA,1),T+1,narray);
idx_sw_cub=zeros(size(FSTA,1),T+1,narray);
eerr_cub=zeros(size(FSTA,1),T+1,narray);


%### UPDATED on 07/14/2020: create variables for market (un)observables
%textFileName_mktvar=sprintf('mktvar.raw');  %the input file 
textFileName_mktvar=sprintf('mktvar_pophhi.raw');  %the input file
mktvar=dlmread(textFileName_mktvar); %the matrix: ahaid year hhi mktcat totpop65
mktvar=sortrows(mktvar,[1,2]);
mktcat_tmp=zeros(size(fs0,1),1);
hhi_tmp=zeros(size(fs0,1),1);
pop65_tmp=zeros(size(fs0,1),1);
for i=1:size(mktvar,1)
    hhi_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,3);
    mktcat_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,4);
    pop65_tmp(fs0(:,1)==mktvar(i,1)&fs0(:,2)==mktvar(i,2),:)=mktvar(i,5);
end
hhi_tmp(fs0(:,3)==0)=0;
mktcat_tmp(fs0(:,3)==0)=0;
pop65_tmp(fs0(:,3)==0)=0;
%hhi=hhi_tmp(:,ones(101,1));
%mktcat=mktcat_tmp(:,ones(101,1));
hhi_cub=repmat(hhi_tmp,1,T+1,narray);
mktcat_cub=repmat(mktcat_tmp,1,T+1,narray);
pop65_cub=repmat(pop65_tmp,1,T+1,narray);


%### UPDATED on 05/05/2020: Add a new input that describes the first next
%three periods' choices, to calculate the 2nd part of sunk cost---the
%ramping part
sunk2=zeros(size(fs0,1),T+1,narray);
DOM=zeros(size(ddom0,1),size(ddom0,2),narray);
eer0=dlmread('eerr1.csv'); %matrix: ahaid year option_given pr1 pr2 ... pr100

%### UPDATED on 05/07/2020: construct the column with the mktdom of the
%current hospital
%{
A=[kron(nowDOM(:,1:2),ones(numel(choice),1)),repmat(choice',size(nowDOM,1),1),reshape(nowDOM(:,3:end)',[],1)];  %matrix A: [ahaid year vid mktdom]
%{
actprob=zeros(size(eer0,1),1);
for k=1:size(eer0,1)
    if eer0(k,3)~=0
       actprob(k)=A(A(:,1)==eer0(k,1)&A(:,2)==eer0(k,2)&A(:,3)==eer0(k,3),end);
    end
end
MKTDOM=[repmat(actprob,1,tbar-1),zeros(size(ddom0,1),T-tbar+2)];
%}
MKTDOM=sortrows(A,[1,2,3]);
MKTDOM(ddom0(:,3)==0,3)=0;
MKTDOM_cub=repmat(MKTDOM(:,end),1,T+1,narray);
clear A
%}


EERR=zeros(size(eer0,1),size(eer0,2)+1,narray);   %#### UPDATED on 05/05/2017: (Take log!!!) the output for the expected logit error: ahaid, option given, prob@t=1, prob@t=2, ..., prob@t=100
clear fs0 ddom0 eer0 
clear hhi_tmp mktcat_tmp pop65_tmp
for c=1:narray
    s=page(c);
    textFileName1=sprintf('fs%d.csv',s);
    textFileName2=sprintf('mks%d.csv',s);   
    textFileName3=sprintf('eerr%d.csv',s);
    tempEERR=dlmread(textFileName3);
    tempEERR(tempEERR==0)=1;
    EERR(:,:,c)=[tempEERR(:,1:3),zeros(size(tempEERR,1),1),log(tempEERR(:,4:end))]; %#### UPDATED on 03/13/2017:include one zero column before all the non-zero expected errors because the first one flow has no expected error.
    EERR(:,:,c)=sortrows(EERR(:,:,c),[1,2,3]);    
    FSTA(:,:,c)=dlmread(textFileName1);  %the matrix: ahaid year bdtot prechs action_given fs1 fs2...fs100
    DOM(:,:,c)=dlmread(textFileName2);  %the matrix: ahaid year action_given ddom_pre ddom0 ddom1 ddom2...ddom99   
    FSTA(:,:,c)=sortrows(FSTA(:,:,c),[1,2,5]);
    %### UPDATED on 05/05/2020: try decreasing sunk cost; FSTA(:,4:6) is
    %for period 0, 1, and simulated period 1; FSTA(:,5:7) is
    %for 1, simulated period 1, and simulated period 2
    fir=FSTA(:,4:(4+tN-t0),c);
    lat=FSTA(:,5:(5+tN-t0),c);
    diff=(fir~=lat);
    [~,loc]=max(diff,[],2);  %loc shows the FIRST index where fir~=lat EXCEPT rows with ALL "0"
    loc(max(diff,[],2)==0)=0; %Let the index for which the rows have ALL zeros be zero
    loc(loc==0)=12;  %Make the index of all-zero rows be 12 so that "sub2ind" can work 
    onepage_sunk2=zeros(size(sunk2,1),T+1);
    idx=sub2ind(size(onepage_sunk2),(1:size(sunk2,1))',loc);
    onepage_sunk2(idx)=(tbar-loc).*(loc>0).*(loc<tbar);
    sunk2(:,:,c)=onepage_sunk2;   
    DOM(:,:,c)=sortrows(DOM(:,:,c),[1,2,3]);
    
    %### UPDATED on 06/22/2020: try contructing cubic so as to remove the for loops in fminsearch 
    % ev of mktdom
    gg_cub(:,:,c)=DOM(:,4:end,c);
    % cost-saving per bed
    fchs_cub(:,:,c)=FSTA(:,4:end-1,c)+11;
    % sunk costs
    fir=FSTA(:,4:end-1,c);
    lat=FSTA(:,5:end,c);
    ind_cub(:,:,c)=(fir~=lat).*lat;
    % switching cost
    idx_sw_cub(:,:,c)=(fir~=lat).*(fir>1);
    % the expected errors
    eerr_cub(:,:,c)=EERR(:,4:end,c);       
end
clear fir lat

%hbed=FSTA(:,3,1); %bdtot: 25%: 25; 50%: 80; 75%: 193
%bdtype=(hbed<=25&hbed>0)+(hbed<=80&hbed>25)*2+(hbed<=193&hbed>80)*3+(hbed>193)*4;


%### UPDATED on 05/07/2020: move this part to the fron to calculate
%absolute mktdom
%{
textFileName26=sprintf('regvar6_sing.raw');  %the input file  
textFileName27=sprintf('regvar7_sing.raw');  %the input file     
textFileName28=sprintf('regvar8_sing.raw');  %the input file     
textFileName29=sprintf('regvar9_sing.raw');  %the input file    

char6=dlmread(textFileName26); %the matrix of characteristics of stand-alone hospitals: ahaid year bdtot prechs adj_vnd hsa
char7=dlmread(textFileName27);    
char8=dlmread(textFileName28);
char9=dlmread(textFileName29);
%}


%%%%%%%%%%%%%%%%%%%%%%%%%%
%with 26 parameters
%theta0=[-6.752912 -6.311207 -3.715860 -4.035594 -6.650967 -6.566256 -7.581982 -5.019611 -5.121039 -5.779016 -4.164671 -4.502513 0.000369 0.000534 0.000058 -0.000058 0.000528 0.000580 0.000505 0.000105 0.000477 0.000506 0.000357 0.000294 0.062659 -1.8];
%theta0=[-5.869513 -5.330296 -2.948035 -3.377029 -5.751234 -5.651111 -6.713429 -4.106970 -4.216721 -4.870577 -3.188701 -3.664923 0.000197 0.000399 -0.000042 -0.000093 0.000394 0.000460 0.000374 -0.000148 0.000345 0.000378 0.000253 0.000127 0.028334 -1.428617];
 

 %theta0=[-6.0418682	-5.7275919	-4.2186812	-3.8558114	-5.7440653	-5.7927766 ...
 % -6.3105595	-4.9054132	-4.9710925	-5.2413514	-4.2134675	-3.7461579 ...
 % 0.000326417	0.000379318	-7.42E-05	-0.001186604	0.000380924	0.000458992 ...
 % 0.000466934	-0.000582455	0.000338908	0.000397496	0.00019452	-5.66E-05 ... 
 % 0.0131	  -0.018  0.027  	0.00043 	-0.0055  	-0.0073 ...	% UPDATED on 05/05/2020: add VARIOUS decreasing sunk cost (2-7)
 % -0.0093	0.0041 	0.00061	 -0.0053	0.0051	 0.000178 ...   % UPDATED on 05/05/2020: add VARIOUS decreasing sunk cost (8-13)
 % -0.0013 ...   %### UPDATED on 05/06/2020: add interaction of desunk and mktdom
 % 0.21580906	0.12910263	0.083096889	0.013831496	-1.0373199	-1.8934718	-0.98035968	-0.27365732];  %### UPDATED on 05/05/2020: add decreasing sunk cost

 
 theta0=[-5.4828365	-5.3015166	-3.7109648	-4.2481047	-5.2400416	-5.5964036 ...
  -5.5251035	-5.0561345	-4.4459762	-4.6060786	-3.9456009	-3.212796 ...
  0.000359496	0.000467073	-0.000135884	0.00010428	0.000407091	0.000483664 ...
  0.000449537	0.000403087	0.000398288	0.000415923	0.000374886	-7.78E-05 ... 
  0.02140	0.17140	0.13873	0.06973	-0.10627	0.00651 ...	% ### UPDATED on 07/14/2020: coefficients for HHI
-0.04973	0.02553	0.17727	0.11473	0.15127	0.16700 ...	
0.00639	0.03547	0.01493	-0.00753	0.00508	-0.00014 ...	% ### UPDATED on 07/14/2020: coefficients for totpop65
-0.02300	0.00213	0.01807	0.02707	0.02067	0.01953 ...	
0.012 0.024 0.6 ...             % ### UPDATED on 07/14/2020: coefficients for market categories
   0 0  0 0 0 0 ... % ### UPDATED on 11/01/2020: add VARIOUS decreasing sunk cost (2-7)
   0 0 0 0 0 0 ...  % ### UPDATED on 11/01/2020: add VARIOUS decreasing sunk cost (8-13)
  0.063134138	0.13487717	0.07500461	0.095128377	-1.753104	-1.8999228	-1.3209065	-0.67105583];  %### UPDATED on 05/05/2020: add decreasing sunk cost



%{
%theta0=[-5.851851	-5.2944001	-2.8233827	-3.196007	-5.7383223	-5.6147304...
 %   -6.8702779	-4.1413453	-4.1940166	-4.7819223	-3.1666896	-3.6302424 ...
  %  1.86681	3.90219	-1.23522  -2.5163	4.01914	4.53675...
  %  4.38982	 -0.996	3.33921	3.28119	2.41559	1.10417...
  %  5.172  2.7550101	-1.49214];   %rescale to reduce rounding error

%theta0=[-5.91740710 -5.37648950 -2.91890680 -3.27645740 -5.79194900 -5.70896380...
  %  -6.77502670 -4.22487880 -4.26353970 -4.90114720 -3.25408260 -3.70270700...
   % 0.00019496 0.00039429 -0.00011557 -0.00023990 0.00038988 0.00045619...
    %0.00037716 -0.00008556 0.00034018 0.00036918 0.00024553 0.00012336...
   % 0.03299912 -1.39623530];
%}
%theta0=[-4.9:-0.1:-6,0.0008:0.0001:0.0019,0.1];
%options = optimset('MaxIter',1e+10000000000000000000000,'MaxFunEvals',1e+10000000000000000000);
options = optimoptions(@fminunc,'MaxIter',1e+10000000000000000000000,'MaxFunEvals',1e+10000000000000000000,'Algorithm','quasi-newton');


%[theta,fval] = fminsearch(@(theta) logL_ReSTAT_decsunkDOM_simp(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,MKTDOM_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,sunk2,theta),theta0,options);
[theta,fval] = fminunc(@(theta) logL_ReSTAT_VARdecsunk(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,sunk2,theta),theta0,options);


textFileName_theta=sprintf('theta_VARdecsunk_fullmkt_w5%d.csv',start);  
dlmwrite(textFileName_theta,[theta fval],'delimiter', ',', 'precision', 8)


%%%%%%%%%%%%%%%%%%
%Compute the standard error
epsilon=1e-6;  %define the disturbance
pr=prob_ReSTAT_VARdecsunk(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,sunk2,theta); %pr26 is the probability function for 26 parameters.
oldpr=pr;
oldlnp=log(pr); %the probability given the original estimates
oldq=theta; %the original estimates
dldq=[];
for i=1:size(theta,2)   %the following for-loop computes the 1st derivative w.r.t. theta
    theta=oldq;
    theta(i)=oldq(i)+epsilon;
    pr1=prob_ReSTAT_VARdecsunk(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,sunk2,theta);
    newlnp=log(pr1);
    dldq=[dldq (newlnp-oldlnp)/epsilon];
end
%A=dldq'*dldq;
%var_q=diag(inv(dldq'*dldq),0);
se_q=sqrt(diag(inv(dldq'*dldq),0));   %the standard error is the sum of the outer product of first derivative

textFileName_se=sprintf('se_VARdecsunk_fullmkt_w5%d.csv',start);  %the input file  
dlmwrite(textFileName_se,se_q,'delimiter', ',', 'precision', 8)
%}

toc
%diary off

